getDictionary <- function(species="mmusculus_gene_ensembl",microarray="affy_mogene_1_0_st_v1",ids) {
	# Get dictionary from BiomaRt, and also TFs from AnimalTFDB
	library(biomaRt)
	mart <- useEnsembl(biomart="ensembl",dataset=species)
	dic <- getBM(attributes=c("ensembl_gene_id",microarray,"external_gene_name"),filters=microarray,values=ids,mart=mart)
	
	TFs <<- scan("../Mus_musculus_TF_EnsemblID.txt",what="character")
	
	return(dic)
}

getExpressionLevels <- function (AllData,experiment) {
	# Print expression data
	expr_file <- paste0("expr-",experiment,".txt")
	
	eset <- rma(AllData)
	e <- exprs(eset)
	dic <<- getDictionary(ids=rownames(e))
	m <- match(rownames(e),dic$affy)
	e <- cbind(dic[m,"external_gene_name"],e)
	write.table(e,expr_file,sep="\t")
	return(eset)
}

differentialExpression <- function(es=eset,des=design,cont=cont.matrix,pval=0.05,experiment) {
	# Find DEGs using limma, and print to file the genes with adj. p-value < pval
	# Calls interally getDictionary
	DE_file <- paste0("DEG-",experiment,".txt")
	
	fit <- lmFit(es, des)
	eb <<- eBayes(contrasts.fit(fit,cont))
	#dic <<- getDictionary(ids=row.names(ebayes))
	DEtable <- topTable(eb,number=Inf,p.value=pval)
	DEgenes <- row.names(DEtable)
	m <- match(DEgenes, dic$affy)
	DEtable$ensembl_gene_id <- dic[m,"ensembl_gene_id"]
	DEtable$gene_name <- dic[m,"external_gene_name"]
	write.table(DEtable,DE_file,sep="\t")
	return(DEtable)
}

filterDEGsByLFCandTFs <- function(DEtable, lfc, tfs=TFs, experiment) {
	# Default: filter for TFs extracted from AnimalTFDB. If tfs=FALSE it is skipped.
	DE_file <- paste0("DE-",experiment,".txt")
	
	DEtable <- DEtable[abs(DEtable$logFC)>log2(lfc),]
	if (class(tfs)=="character") {
		DEtable <- DEtable[which(DEtable$ensembl_gene_id %in% tfs),]
		}
	write.table(DEtable,DE_file,sep="\t")
	return(DEtable)
}

makePhenoMatrix <- function(DEtable,experiment) {
	# Booleanize expression of DEGs
	pheno_file <- paste0("pheno-",experiment,".txt")
	
	vars <- eb$s2.post[rownames(DEtable)] #add variance as calculated by ebayes function
	v <- cbind(DEtable,vars)
	v <- v[order(v$gene_name,-v$vars),]
	v <- v[! duplicated(v[,"gene_name"]),] #keep only the probe with higher variance
	
	pheno_matrix <- cbind(v$gene_name,as.integer(v$logFC >0))
	pheno_matrix <- pheno_matrix[! is.na(pheno_matrix[,1]),] #remove probes which didn't map to any gene
	write.table(pheno_matrix,pheno_file,sep="\t",row.names=F,col.names=F,quote=F)
	return(pheno_matrix)
}

fromHairballToAdjMatrix <- function(pheno_matrix, syn=synonyms, int=interactions, experiment) {
	# Takes MetaCore files and returns the adjacency and phenotypic matrices
	adj_matrix_file=paste0("adjmat-",experiment,".txt")
	red_pheno_file=paste0("phenored-",experiment,".txt")
	
	if (all(rownames(pheno_matrix)==1:nrow(pheno_matrix))){
		rnames=pheno_matrix[,1] 
	} else {
		rnames=rownames(pheno_matrix)
	}
	
	adj_matrix <- matrix(0,nrow=nrow(pheno_matrix),ncol=nrow(pheno_matrix),dimnames=list(rnames,rnames))	#simmetric matrix of 0s
	for (i in (1:nrow(int))) {
		input <- syn[syn[,1]==int[i,1],2]	#this may result in multiple inputs
		output <- syn[syn[,1]==int[i,3],2]	#this may result in multiple outputs
		value <- switch(int[i,2],"-->"=1,"--|"=-1,"--?"=2)
		for (j in input){	#I will assign the same interaction to all the genes corresponding to the network objects...
			for (k in output){
				if ((j %in% rownames(adj_matrix)) && (k %in% rownames(adj_matrix))) { #...if they are DE
					adj_matrix[j,k] <- value
				}
			}
		}
	}
	zeros <- intersect(which(rowSums(abs(adj_matrix))==0),which(colSums(abs(adj_matrix))==0)) #genes without in or out-going interactions
	adj_matrix_red <- adj_matrix[-zeros,-zeros]
	#write.table(adj_matrix_red,adj_matrix_file,sep="\t",quote=F)
	
	pheno_matrix_red <<- pheno_matrix[which(rnames %in% rownames(adj_matrix_red)),]
	#write.table(pheno_matrix_red, red_pheno_file, sep="\t",row.names=TRUE,col.names=F,quote=F)
	
	return(adj_matrix_red)
}

adjMatToTranslatedInteractions = function(adj_matrix,experiment) {
	# Rewrites MetaCore interactions to be comparable with contextualized ones
	trans_int_file=paste0("tint-",experiment,".txt")
	n=rownames(adj_matrix)
	x <- c(); y <- c(); z <- c();
	for (i in (1:length(n))){
		for (j in (1:length(n))){
			if (adj_matrix[i,j] !=0) {
				value <- switch(as.character(adj_matrix[i,j]),"1"="-->","-1"="--|","2"="--?")
				x <- c(x,n[i])
				y <- c(y,value)
				z <- c(z,n[j])
			}
		}
	}
	write.table(data.frame(x,y,z),trans_int_file,row.names=F,col.names=F,quote=F)
}

### Satoshi's variance
select_unique_probe = function(M,names){
	variance = apply(M,1,var,na.rm=TRUE) 
	names(variance)=names
	variance.aggregated.index =tapply(seq_along(variance), names(variance), function(x){
		m = max(variance[x])
		I = which(variance[x]==m)
		x[I[1]]
		})
	data = M[variance.aggregated.index,]
	rownames(data) = names[variance.aggregated.index]
	data
}

# tint->adjmat
tintToAdjmat=function(tint_file,adjmat_file){
	grn=read.table(tint_file,sep=" ",quote="",stringsAsFactors=FALSE)
	colnames(grn)=c("FROM","Effect","TO")

	#Satoshi's script
	enum = matrix(c(seq_along(unique(c(grn[,"FROM"],grn[,"TO"]))), unique(c(grn[,"FROM"],grn[,"TO"]))), ncol=2)
	b = cbind(enum[match(grn[,"FROM"],enum[,2]),1],enum[match(grn[,"TO"],enum[,2]),1]); class(b)='numeric'
	length(unique(c(b[,1],b[,2]))) #unique node number 
	b = cbind(b,grn[,c("FROM","TO")])
	num = length(unique(c(as.numeric(b[,1]), as.numeric(b[,2]))))
	mat = matrix(0,num,num)
	rownames(mat)=colnames(mat)=enum[,2]
	# 1. Indices of activation
	act = grn[grn[,"Effect"]=='-->',]
	mat[cbind(match(act[,1],rownames(mat)),match(act[,3],colnames(mat)))]=1
	# 2. Indices of inhibition
	inh = grn[grn[,"Effect"]=='--|',]
	mat[cbind(match(inh[,1],rownames(mat)),match(inh[,3],colnames(mat)))]=-1
	# 3. Indices of unspecified
	uns = grn[tolower(grn[,"Effect"])=='--?',]
	mat[cbind(match(uns[,1],rownames(mat)),match(uns[,3],colnames(mat)))]=2

	zeros <- intersect(which(rowSums(abs(mat))==0),which(colSums(abs(mat))==0))
	if (length(zeros)>0) mat <- mat[-zeros,-zeros]
	write.table(mat,adjmat_file,sep="\t",quote=FALSE)
	return(mat)
}

#  tint->adjmat, without files
simpleTintToAdjmat=function(grn){
	colnames(grn)=c("FROM","Effect","TO")
	#Satoshi's script
	enum = matrix(c(seq_along(unique(c(grn[,"FROM"],grn[,"TO"]))), unique(c(grn[,"FROM"],grn[,"TO"]))), ncol=2)
	b = cbind(enum[match(grn[,"FROM"],enum[,2]),1],enum[match(grn[,"TO"],enum[,2]),1]); class(b)='numeric'
	length(unique(c(b[,1],b[,2]))) #unique node number 
	b = cbind(b,grn[,c("FROM","TO")])
	num = length(unique(c(as.numeric(b[,1]), as.numeric(b[,2]))))
	mat = matrix(0,num,num)
	rownames(mat)=colnames(mat)=enum[,2]
	# 1. Indices of activation
	act = grn[grn[,"Effect"]=='-->',]
	mat[cbind(match(act[,1],rownames(mat)),match(act[,3],colnames(mat)))]=1
	# 2. Indices of inhibition
	inh = grn[grn[,"Effect"]=='--|',]
	mat[cbind(match(inh[,1],rownames(mat)),match(inh[,3],colnames(mat)))]=-1
	return(mat)
}